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with impact velocity dependent restitution coefficient 
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We report numerical simulations of strongly vibrated granular materials designed to mimic recent 
experiments performed both in presence or absence of gravity. The coefficient of restitution used here 
depends on the impact velocity by taking into account both the viscoelastic and plastic deformations 
of particles, occurring at low and high velocities respectively. We show that this model with impact 
velocity dependent restitution coefficient reproduce results that agree with experiments. We measure 
the scaling exponents of the granular temperature, collision frequency, impulse, and pressure with 
the vibrating piston velocity as the particle number increases. As the system changes from a 
homogeneous gas state at low density to a clustered state at high density, these exponents are all 
found to decrease continuously with the particle number. All these results differ significantly from 
classical inelastic hard sphere kinetic theory and previous simulations, both based on a constant 
restitution coefficient. 

PACS numbers: 05.45.Jn, 05.20.Dd, 45.70.-n 
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I. INTRODUCTION 

The past decade has seen the publication of many ex- 
perimental 0. ■ numerical 000. an( i theoretical 
EI 13 studies of strongly vibrated granular media. 
This problem is interesting because vibrated granular me- 
dia are simple but nontrivial example of a noncquilibrium 
steady states and the only way to experimentally realize 
granular gases Q. However, numerous questions remain 
about the link between experiments on one hand, and 
theory and simulations on the other. Most numerical 
and theoretical studies were not intended to be compared 
with experiments. Therefore, they have parameter values 
far from the experimental ones, and none of them predict 
even the most basic features of the experimental results. 

In this paper, we bridge the gap between experiments 
and numerics by presenting simulations of strongly vi- 
brated granular materials designed to mimic recent ex- 
periments performed both in presence |10| or absence pd| 
of gravity. We present the first simulations which resem- 
ble the experiments for a large range of parameters. We 
show that two parameters are especially important for 
the agreement between experiment and simulation. First 
of all, the coefficient of restitution has to be dependent 
on the particle impact velocity by taking into account 
both the viscoelastic and plastic deformations of particles 
occurring at low and high velocities respectively. Most 
previous numerical studies consider only constant resti- 
tution coefficient 0,0,0], or few studies with the slightly 
velocity dependence (due the only viscoelastic contribu- 
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tion) |l2j . Secondly, it is important to explicitly consider 
the number of particles N. Studying only one value of 
N or comparing results obtained at different iV can lead 
to interpretive difficulties. 

Beyond these agreements between experiments and our 
simulations, we find new results that differ significantly 
from classical inelastic hard sphere kinetic theory and 
previous simulations. We measure the scaling exponents 
of the granular temperature, collision frequency, impulse, 
and pressure with the vibrating piston velocity as the par- 
ticle number increases, both in the presence or absence 
of gravity. We show that the system undergoes a smooth 
transition from a homogeneous gas state at low density 
to a clustered state at high density. 

The paper has the following structure. In Sec. [HJ 
we present a description of the simulations (notably the 
model of impact velocity dependent restitution coeffi- 
cient, and the influence of other simulations parameters). 
Section IlIII provides a comparison of simulations and ex- 
periments (showing the importance of the variable coef- 
ficient of restitution and the particle number), and the 
results of the scaling exponents. Section fill CI focus on 
the influence of other simulations parameters (bed height, 
box size, particle rotations, gravity). Finally, in Sec. II VI 
we summarize our results. 



II. DESCRIPTION OF THE SIMULATIONS 

A. The variable coefficient of restitution 

The greatest difference between our simulations and 
the previous numerical studies of vibrated granular me- 
dia 0, is that we use a restitution coefficient that 
depends on impact velocity. The restitution coefficient r 



Typeset by REVTeX 



2 



is the ratio between the relative normal velocities be- 
fore and after impact. In all previous simulations of 
strongly vibrated granular media, the coefficient of resti- 
tution is considered to be constant and lower than 1. 
However, since a century, it has been shown from im- 
pact experiments that r is a function of the impact ve- 
locity v Q El E E3- Indeed, for metallic particles, 
when v is large (v > 5 m/s 0]), the collidin g p arti- 
cles deform fully plastically and r oc d -1 / 4 [l3l IT3. HI}. 
When u < 0.1 m/s 0], the deformations are elastic 
with mainly yiscoelastic dissipation, and 1 — r oc v 1 / 5 
[TH Il6l UtI Il8| . Such velocity-dependent restitution coef- 
ficient models have rece ntly shown to be important in nu- 
merical [HE |IJ,|2l|, HIE! and experimental [HIH 
studies. Applications include: granular fluidlike proper- 
ties (convection 0], surface waves pojl, collective col- 
lisional processes (energy transmissio n l2l| , absence of 
collapse E,H2) and planetary rings |23l l24j . But sur- 
prisingly, such model has not yet been tested numerically 
for strongly vibrated granular media. 

In this paper, we use a velocity dependent restitution 
coefficient r(v) and join the two regimes of dissipation 
(viscoelastic and plastic) together as simply as possible, 
assuming that 




(m/s) 



FIG. 1: The restitution coefficient r as a function of impact 
velocity v, as given in Eq. (0 (solid line). The dashed lines 
show vq = 0.3 m/s and ro = 0.95. Experimental poin ts (•) 
for steel spheres were extracted from Fig.l of Ref. 25] 
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v < v , 
V > v , 



(1) 



where vq = 0.3 m/s is chosen, throughout the paper, to 
be the yielding velocity for stainless steel particles |14 . l25| 
for which r is close to 0.95 p5|. Note that vq ~ 1/y/p 
where p is the density of the sphere 0|. We display 
in Fig. ^ the velocity dependent restitution coefficient of 
Eq. 0, with r = 0.95 and vo — 0.3 m/s, that agrees 
well with experimental results on steel spheres from Ref. 
[25). As also already noted by Ref. Q, the impact ve- 
locity to cause yield in metal surfaces is indeed relatively 
small. For metal, it mainly comes from the low yield 
stress value (Y ~ 10 9 N/m 2 ) with repect to the elastic 
the Young modulus (E ~ 10 11 N/m 2 ). Most impacts 
between metallic bodies thus involve some plastic defor- 
mation. 



B. The other simulation parameters 

The numerical simulation consists of an ensemble of 
identical hard disks of mass m«3xl0~ 5 kg excited ver- 
tically by a piston in a two-dimensional box. Simulations 
are done both in the presence (g — 9.8 m/s 2 ) and absence 
(g = 0) of uniform gravity g. Collisions are assumed in- 
stantaneous and thus only binary collisions occur. For 
simplicity, we neglect the rotational degree of freedom. 
Collisions with the wall are treated in the same way as 
collisions between particles, except the wall has infinite 



Motivated by recent 3D experiments on stainless steel 
spheres, 2 mm in diameter, fluidized by a vibrating pis- 
ton |Toj. we choose the simulation parameters to match 
the experimental ones: in the simulations, the vibrated 
piston at the bottom of the box has amplitude A = 25 
mm (distance between the highest and lowest positions 
of the piston) and frequencies 5 Hz < / < 50 Hz. The 
piston is nearly sinusoidally vibrated with a waveform 
made by joining two parabolas together. The vertical 
displacement of the piston z(t) during time t then is 
z(t) = 4r{t 2 -tl) for -t < t < t and z(t) = -^(t 2 -t%) 
for to < t < 3to with t = 1/(4/). This leads to a max- 
imum piston velocity given by V = AAf '. The particles 
are disks d = 2 mm in diameter with stainless steel colli- 
sion properties through v and r (see Fig. 0. The box 
has width L = 20 cm and horizontal periodic boundary 
conditions. Since our simulations are two dimensional, 
we consider the simulation geometrically equivalent to 
the experiment when their number of layers of particles, 
n = Nd/L, are equal. Hence in the simulation, a layer 
of particles, n = 1, corresponds to 100 particles. We 
checked that n is an appropriate way to measure the num- 
ber of particles by also running simulations at L = 10 cm 
and L = 40 cm. None of this paper's results depend sig- 
nificantly on L. As in the experiments, the height h of the 
box depends on the number of particles in order to have 
a constant difference h — ho = 15 mm, where ho is the 
height of the bed of particles at rest. Heights are defined 
from the piston at its highest position. The influence of 
h— ho on the results is discussed in Sec. IIII CI 



3 



(a) experiments (b) r = r(v) 




(c) r = 0.95 (d) r = 0.7 




FIG. 2: Time averaged pressure P on the top of the cell as 
a function of particle layer, n, for various vibration frequen- 
cies, / : (a) Experimental results from [To| for stainless steel 
beads 2 mm in diameter, with A = 25 mm, 10 < / < 20Hz 
with a 1 Hz step (from lower to upper) and h — ho = 5 mm. 
(b) Numerical simulation where the coefficient of restitution 
is given by Eq. (c) Numerical simulation with a coef- 
ficient of restitution is 0.95, independent of impact velocity, 
(d) Numerical simulation with a coefficient of restitution is 
0.7, independent of impact velocity. The simulations (b) - 
(d) are 2D with gravity, done for 2 mm disks, with A — 25 
mm, 10 < / < 30 Hz with a 2 Hz step (from lower to upper) 
and h — ho = 15 mm. In the simulations, the two-dimensional 
pressure is given in N/m. 



III. COMPARISON OF SIMULATION AND 
EXPERIMENT - SCALING PROPERTIES 

A. The importance of the variable coefficient of 
restitution 

We examine first the dependence of the pressure on 
the number of particle layers for maximum velocity of 
the piston 1 < V < 5 m/s (V = AAf). The time av- 
eraged pressure at the upper wall is displayed in Fig. |3 
as a function of n for various /: from the experiments 
of Falcon et al. [13 (see Fig. [It), fro m our simulations 
with velocity dependent restitution coefficient r = r(v) 
proposed in Eq. (see Fig.[2b), and with constant resti- 
tution coefficient r = 0.95, often used to describe steel 
particles (see Fig. [2:), and finally with an unrealistic con- 
stant restitution coefficient r = 0.7 (see Fig. |2Ji). Simu- 
lations with r = r(v) give results in agreement with the 
experiments: At constant external driving, the pressure 
in both Figs. [3]} and|2jD passes through a maximum for a 
critical value of n roughly corresponding to one particle 
layer. For n < 1, most particles are in vertical ballistic 
motion between the piston and the lid. Thus, the mean 



pressure increases roughly proportionally to n. When n 
is increased such that n > 1, interparticle collisions be- 
come more frequent. The energy dissipation is increased 
and thus the pressure decreases. This maximum pressure 
is not due to gravity because it also appears in simula- 
tions with ,g = and r = r(v). Furthermore, the max- 
imum persists when g is increased above 9.8 m/s 2 . For 
n > 4 and for certain frequencies, a resonance appears in 
Fig. I2t> which is controlled by the ratio between the vi- 
bration period and the particle flight time under gravity, 
\fgfhj f ■ Turning our attention to Fig. we see that 
setting r — 0.95 independently of impact velocity gives 
pressure qualitatively different from experiments. The 
difference between Figs. 03 andEt can be understood by 
considering a high velocity collision (e.g. v = 1 m/s). 
In Fig. [2b, this collision has a restitution coefficient of 
r = r(l m/s) w 0.7 (see Fig. [I), whereas in Fig. |2J;, r 
is fixed at 0.95 for all collisions. This means that for 
equal collision frequencies, dissipation is much stronger 
for r = r(v) than for r = 0.95, because the high velocity 
collisions dominate the dissipation. Stronger dissipation 
leads to lower granular temperatures and thus to lower 
pressures. 

We can check this interpretation by changing the con- 
stant restitution coefficient to r = 0.7 and then compar- 
ing it to r — r(v). In these two cases, the high velocity 
collisions will have roughly the same restitution coeffi- 
cient. We indeed observed a pressure that decreases for 
large n for constant r — 0.7 (see Fig. |2Ji). Therefore, 
surprisingly, constant r = 0.7 reproduces more precisely 
the experimental pressure measurements than constant 
r = 0.95, even though r — 0.95 or r = 0.9 is often given 
as the restitution coefficient of steel. However, if we look 
at other properties, we see that r = 0.7 and r = r(v) give 
very different predictions. 



(a) r — r(v) 




FIG. 3: Snapshots from the simulations with n = 3, gravity 
j / 0, driving frequency / = 30 Hz and h — ho = 15 mm. 
The upper wall is stationary, and the lower wall is the piston, 
which is at its lowest position in both snapshots. The hor- 
izontal boundaries are periodic (indicated by dashed lines). 
Gravity points downwards, (a) r — r(v), as given in Eq. Q, 
and (b) constant r = 0.7. In (b) we see a tight cluster which 
was not observed in the experiments. 



4 



For example, in Fig. |3 we show two snapshots from 
two different simulations one with r — r(v) and another 
with r = 0.7, both with n = 3 in the presence of grav- 
ity. When r = r(v), the particles are concentrated in the 
upper half of the chamber, but they are evenly spread 
in the horizontal direction (see Fig. |3K) - The system is 
hotter and less dense near the vibrating wall, and colder 
and denser by the opposite wall. But, when r = 0.7, the 
majority of the particles are confined to a tight cluster, 
pressed against the upper wall, coexisting with low den- 
sity region (see Fig-Efc) ■ This instability has been already 
been reported numerically |2t| , although for much differ- 
ent parameters (constant restitution coefficient r = 0.96, 
thermal walls, no gravity, and large n). However, nothing 
like this was ever seen experimentally. Therefore, if one 
is seeking information about particle positions, r = 0.7 
gives incorrect results even though it gives acceptable re- 
sults for the pressure. We conclude, therefore, that the 
only way to successfully describe all the properties in all 
situations is to use a velocity dependent restitution coef- 
ficient model. 



B. The importance of the particle number 

Many authors have postulated that the pressure on the 
upper wall P (or granular temperature T) is related to 
the piston velocity V through P, T cx V . However, it is 
not clear what the correct "scaling exponent" 8 should 
be. This question has been addressed several times in the 
past, without a clear resolution of the question HUES HI- 
For example, kinetic theory P, and hydrodynamic 
models jgl predict T cx V 2 whereas numerical simula- 
tions 0, 3 or experiments 0, Q gi ye T cx V 6 , with 
1 < 8 < 2. These studies were done at single values of 
n. In this section, we show that it is very important to 
explicitly consider the dependence of the scaling expo- 
nents on n. We also consider the effect of gravity and a 
variable coefficient of restitution. Doing so enables us to 
explain and unify all previous works. 

At the upper wall, we measured numerically the col- 
lision frequency, N c and the mean impulsion per colli- 
sion, AI for various frequencies of the vibrating wall and 
numbers of particles in the box, with r = r(v) or with 
r = 0.95, in the presence or absence of gravity. The time 
averaged pressure on the upper wall can be calculated 
from these quantities using 

P = N C AI/L . (2) 

(By conservation of momentum, the time averaged pres- 
sure on the lower wall is just P plus the weight of the 
particles Nmg/L.) The total kinetic energy of the sys- 
tem is also measured to have access to the granular tem- 
perature, T. N c , AI, P, and T are all found to fit with 
power laws in V 8 for our range of piston velocities. Fig. 
0] shows 8 exponents of N c , AI, P, and T as a function 
of n. When g = and r is constant (see Fig. 2b) > we 
have P ~ V 2 , AI ~ V and N c ~ V for all n. We call 



these relations the classical kinetic theory scaling. This 
scaling can be established by simple dimensional analy- 
sis when the vibration velocity V provides the only time 
scale in the system. This is the case for g = and r in- 
dependent of velocity. However, in the experiments, two 
additional time scales are provided, one by gravity and 
another by velocity dependent restitution coefficient. Nu- 
merical simulations can separate the effects of these two 
new time scales on the scaling exponents 8. This is done 
in Fig. 2t> [where g = but r = r(v)] and Fig. 2b [where r 
is constant but g ^ 0] . In both figures, all the exponents 
become functions of n. However, the time scale linked 
to r = r(v) leads to much more dramatic departure from 
the classical scaling. After considering the two time scale 
separated, let us consider the case corresponding to most 
experiments, where both gravity and r = r(v) are present 
(see Fig. The similarity between this figure and Fig. 
2b confirms that the velocity dependent restitution co- 
efficient has a more important effect than the gravity. 
Furthermore, only the variation of the restitution coef- 
ficient on the particle velocity explains the experiment 
performed in low-gravity This experiment gives a 

y3/2 p ressure scaling (•-mark on Fig. 2b) for n — 1 and 
a motionless clustered state for n > 2. Only the simula- 
tion with r = r(y) can reproduce these results (see Fig. 
0b) whereas constant r simulations leads to the classical 
scaling (P cx V 2 , see Fig. 2b) an d only a gaseous state 
for all n shown in the figure. 



As shown in Fig 21 it is thus very important to ex- 
plicitly consider the dependence of 8 on n. In all cases, 
except the unrealistic case of Fig. 2bj depends on n. 
To our knowledge, the only experiment (To) to systemat- 
ically investigate this effect shows that T cx V etyn \ with 
8 continuously varying from 8 — 2 when n — > 0, as ex- 
pected from kinetic theory, to 8 ~ for large n due to 
the clustering instability. These experiments |l0j | per- 
formed under gravity (shown in Fig. [SJ are well repro- 
duced by the simulations of Fig. 2b- I n both cases, the 
observed pressure and granular temperature scaling ex- 
ponents strongly decrease with n. 



We finish this section by noting two curious facts about 
Fig. 21 First of all, in Fig. 2b [g = and r = r(v)}, 9ml 
for the pressure and temperature when n > 2. This is 
the sign of new robust scaling regime where P and T cx 
V , which will be the topic of a future paper. Secondly, 
in Fig. 2b \ff 7^ and r = r(v)], the scaling exponents 
are not shown for n > 3, because the dependence of P, 
T, N c and AI on V is no longer a simple power law. 
(More precisely, we do not plot a point on Fig. 21 when 
|logX obs0 rvcd - logAfittcdl > 0.25 for any of the eleven 
simulations used to calculate the exponent - see caption.) 
The power law breaks down because there is a resonance 
between the time of flight of the cluster under gravity 
and the vibration period. 
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(a) r = 0.95, g = 

2| ge0GO0G00Q Q Q ' Q O p 0-0*0 QQO O O O p O O O O p 



(b) r = r(«), 3 = 




FIG. 4: The exponents # as a function of n which give the scaling of the granular temperature T (0), collision frequency N c 
(*), mean impulsions AI (A) and pressure P (o). All these quantities are proportional to V 6 ^ n K Without gravity: (a) for 
r — 0.95 and (b) for r = r(v). With gravity: (c) for r = 0.95 and (d) for r = r(v). The exponents are obtained by fixing n 
and performing eleven simulations, varying / from 10 Hz to 30 Hz. Then logX (where X is the quantity being considered) is 
plotted against log V. The resulting curve is always nearly a straight line (except for n > 3 in (d) - see text), and the exponent 
is calculated from a least squares fit. The pressure scaling point (•) on (b) is from experiment |Tl| performed in low-gravity. 
See Fig. (resp. Fig.[3p) for typical snapshots corresponding to n = 3, g ^ and r = r(v) (resp. r = 0.95) 
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FIG. 5: Experimental data performed under gravity from Ref. 
[l(if : The exponents #(n) of time averaged pressure (□) (see 
Fig-Et)i an d kinetic energy extracted from density profile (o) 
or volume expansion (0) measurements. These data should 
be compared with the simulations of Fig. 0)1. 



C. The influence of other parameters 

In this section, we review the influence of the other 
simulation parameters (box size, particle rotations, grav- 



ity, and the vibration parameters), and show that it is not 
possible to reproduce the experimental curves in Fig. [5^, 
unless one sets r — r(v) or r = 0.7. 

Performing simulations for 5 mm < h — ho < 50 mm 
show that the shapes of the curves P vs. n in Fig. [2b and 
[3: remain the same. For r — r(v) (Fig. increasing 
h—h,Q shifts the maximum towards smaller values of n and 
decreases in amplitude. The only exception occurs when 
box height approaches the particle diameter, i.e. h — ho = 
5 mm, where the maximum disappears. Considering r — 
0.95 leads to similar conclusions. 

To eliminate the possibility that experimental curve 
can be reproduced by taking into account particle ro- 
tations, we performed simulations with r = 0.95 and 
various values of the tangential restitution coefficient r t . 
This parameter is defined as the ratio between the tan- 
gential components of the pre- and post-collision rela- 
tive velocities. Perfectly smooth spheres correspond to 
r t = — 1. When r t = 1, the tangential relative velocity is 
reversed by the collision. These two values, \r t \ = 1, cor- 
respond to energy conservation. Energy is dissipated for 
— l<7"t<l, rt = corresponding to maximum energy 
dissipation. When \r t \ is close to one, the P vs. n curves 
are almost unchanged. When r t is close to 0, the curves 
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become nearly flat for n > 2. 

Throughout this paper, we have used the piston vi- 
bration velocity V to characterize the vibration, ft is 
important to point out that V is not the only way to do 
this. One could also use the maximum piston accelera- 
tion r. When r is close to g, it controls the behavior of 
the system, i.e., adjusting A and / while keeping T con- 
stant does not change the system's behavior much. But 
in the simulations presented here, r> j, and the sys- 
tem's behavior is controlled by V. This can be checked 
by multiplying the frequency by 10 while dividing A by 
10, thus keeping V the same (while T increases by an 
order of magnitude) . Doing so changes the pressure only 
by about 20%. Therefore, V is the correct parameter 
to describe the vibration for the simulations considered 
here. 



IV. CONCLUSIONS 

In this paper, we brought simulations of a strongly 
vibrated granular medium as close as possible to the ex- 
periments. We show that the use of a velocity depen- 
dent coefficient of restitution reproduce results that agree 
with experiments. It is especially important to take into 
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